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The initial stages of relativistic heavy ion collisions are studied numerically in the framework of 
a 2+1 dimensional classical Yang-Mills theory. We calculate the energy and number densities and 
momentum spectra of the produced gluons. The model is also applied to noncentral collisions. The 
numerical results are discussed in the light of RHIC measurements of energy and multiplicity and 
other theoretical calculations. Some problems of the present approach are pointed out. 
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I. INTRODUCTION 

In ultrarelativistic heavy ion collisions such as stud- 
ied at RHIC and LHC particle production in the central 
rapidity region is dominated by the gluonic degrees of 
freedom in the nucleus. At sufficiently small x the phase 
space density of these gluons is large, so one can try to 
treat them as a classical color field. Let us first briefly 
review the model of 0, 0, before turning to our results 
in Sec. I1VI and their phenomenological implications in 
Sec. E] Our notation is essentially that of 0, El • 

The idea of was to model the high momentum de- 
grees of freedom of a nucleus as static random classical 
color sources with a Gaussian probability distribution: 

(p a (x T )p b (y T )) = g 2 p 2 S ab S 2 (x T - y T ), (1) 

where xt and yx are vectors in the transverse plane. 
The classical color field generated by this source is then 
obtained from the equations of motion 

[D^F^] = J V . (2) 

This original formulation of the model is very simple, 
beyond the nuclear radius Ra it only depends on one di- 
mensionful phenomenological parameter p, (related to A s 
introduced in [4( by A s = g 2 p) and the QCD coupling g 
that does not run in this classical approximation. One 
may, however, argue that the Gaussian probability dis- 
tribution should be replaced by something else, namely a 
solution of the "JIMWLK" renormalization group equa- 
tion @. 

The McLerran-Venugopalan model [jj describes the 
wave function of one nucleus. Nucleus-nucleus collisions 
were first studied in this framework in Q. The source 
current is taken to be 

J" = 5» + p (1) (x T )5(x-) + 6»-p (2) (x T )6(x+), (3) 

where the color charge densities P{ m ) of the two nuclei 
are independent. In the region x~ < 0, x + < which is 
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causally connected to neither of the nuclei, the solution 
can be chosen as = 0. In the regions x~ < 0, x + > 
and x~ > 0, x + < which are causally connected to 
only one of the nuclei the solutions are "transverse pure 
gauges" 

A\ m) = _V A <">cfe- jA <->, with V|A (m) = -gp (m) . 

(4) 

The initial condition (r = 0) for the interesting region 
x~ > 0, x + > is obtained by matching the solutions on 
the light cone. This yields: 

A'| T=0 = A\ 1)+ A\ 2) , (5) 

■^■ ?? |r=Q = Y^C 1 )' ^(2)1- 

Modeling the sources as delta functions on the light 
cone (Eq. J3J|) makes the initial conditions boost invari- 
ant. We shall also restrict ourselves to strictly boost in- 
variant field configurations. This elimination of the longi- 
tudinal degrees of freedom makes the numerical solution 
of the equations of motion easier, but is a serious lim- 
itation, especially for studying thermalisation (see e.g. 
0). 



II. (2+l)-DIMENSIONAL CLASSICAL 
HAMILTONIAN CHROMODYNAMICS ON THE 
LATTICE 

The analytic solution of the equations of motion, Eq. 
@, with the initial conditions, Eqs. J5J, is not known, 
but they can be studied numerically. A lattice Hamil- 
tonian formulation of the model was first developed in 
0. 

Assuming that the field configurations are boost- 
invariant reduces the system to a 2+1-dimensional one. 
Choosing the Schwinger gauge A T — one can cast the 
equations of motion into a Hamiltonian form. The lattice 
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Hamiltonian is 

= E 



9 

y a 
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- — Tr 7r 2 
a 



a \ -> „ 



(6) 



where a is the lattice spacing and Ei,Ui,ir and <f> are 
dimensionless lattice fields. The fields are matrices in 
color space, with E % = E a l t a etc. and the generators 
of the fundamental representation normalised in the con- 
ventional way as Tr t a tb = l/2S a b- The first two terms 
are the transverse electric and magnetic fields, with the 
transverse plaquette 

U ± {x T ) = U x (x T )U y (x T + e x )U x (x T + e y )Ul{x T ). (7) 

The last two terms are the kinetic energy and covariant 
derivative of the rapidity component of the gauge field 
(f> = A, q = — t 2 A ?? , which becomes an adjoint represen- 
tation scalar with the assumption of boost invariance. 
For the parallel transported scalar field we have used the 
notation 



(j>i(x T ) = Ui(x T )4>(xT + ei)Uj(x T ) 



(8) 



In the Hamiltonian, Eq. 10, there is a residual invari- 
ance under gauge transformations depending only on the 



transverse coordinates. The Hamiltonian equations of 
motion arc 



Us = 



E x = 



E y 
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T 
T7T 
IT 

i rT 



E l Ui (no sum over i), 



7^2 \Ux, y + U x - y - h.c. ] - trace 



T 
IT 
2 



(9) 
(10) 

(11) 



2.9 



h.c. 



trace 



+ 4>-i — 20 



(12) 



The Gauss law, conserved by the equations of motion, 
reads: 

\u\(x T - ei)E\x T - ei)Ui{x T - e*) - E\x T ) 
- #,tt]=0. (13) 



On the lattice the initial conditions (0) become 

Tr 





E' 

■k{x t ) 



t a ((uP + u? : 



1 + Uj)- h 



^£[(^)-l) (ul^{x T )-U^(x T )) 
{u}{x T - a) - l) (uf >(s T - ei) - UP (xr ~ *)) 



(14) 

(15) 
(16) 

(17) 



- h. 



where L^ 1 ' 2 ) in Eq. (|14(l are the link matrices correspond- 
ing to the color fields of the two nuclei (A[ 1 ' 2 ' ) in Eq. Q) 
and the link matrix Ui corresponding to the r > color 
field Ai must be solved from Eq. (|14|) . 

The model has three free parameters, the coupling g, 
the source density [i and the nuclear transverse area ttR\. 
In this work the lattice size is taken to be L 2 — N 2 a 2 — 
ttR 2 a . This means that the field modes have an infrared 
cutoff of the order 1/Ra, while physically one would ex- 
pect them to be cut off at a scale ~ Aqcd by confinement 
physics not included in the classical field model. So in 
order to be physically sensible our results should not de- 
pend on this infrared cutoff. 



The values of the three parameters g, [i and kR 2 a 
separately are needed when translating lattice units to 
physical units, but the dimensionless parameter g 2 [iRa 
controls the qualitative behavior of the model; the weak 
coupling or weak field limit is reached for small values 
of this parameter (see also Q). To see this consider the 
system on a transverse lattice of spacing a. Now we have 
S 2 (xt) ~ l/a 2 . Thus, from Eq. (JIJ, the charge den- 
sity p ~ ,9/i/a. The Green's function of the operator V 2 ^ 
in Eq. is a logarithm, which is parametrically con- 
stant. Thus A(cct) is obtained by summing contributions 
~ g-a 2 ■ gji/a from each of the ~ R\j a 2 cells (the area of 
a cell being a 2 ). Because the charges are distributed as 
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Gaussians with zero expectation value, their sum scales 
as a square root of the number of lattice sites, and we 
get A(xt) ~ ag 2 p ■ \J R\ / a 2 ~ g 2 \iRA- Because of the 
exponentials of A in Eq. J3J it is the magnitude of the 
dimensionless field A that determines the nonlinearity of 
the model. 

The same argument can also be formulated in mo- 
mentum space. The Poisson equation, Eq. (J3J, can 
be written as kT 2 A(kx) = gpik^)- One needs a pre- 
scription to deal with the zero mode, the one chosen 
here is color neutrality of the system as whole, p{kx = 
Ot) = = A(fcr = Ot)- Then the dominant contri- 
bution comes from the smallest nonzero Fourier mode, 
k T ~ 1/Ra- In momentum space the correlator Q is 
{p(k T )p{p T )) 
Thus A(fe T ) r 



g 2 p?5 2 {k T + Pt) with 5 2 {k T ) 



gR\p(k T ) 



gR A ■ gp.RA and A(xt) 



A(k T )/R 2 A ~ g 2 \xR A - 

In this Hamiltonian formalism the energy per unit 
rapidity in different field components is naturally the 
easiest and the most fundamental quantity to compute. 
One can also measure equal time correlation functions of 
fields: 



(E?(k T ,T)E?(~k T ,T)), 

(A«(k T ,T)AU-k T ,T)), 

(vr a (fc T ,rK(-fcT,T)), 

(r(kT,r)r(-k T ,r)). 



(18) 
(19) 
(20) 
(21) 



These correlation functions are not gauge invariant. One 
can, however, argue that in the Coulomb gauge diAi = 
a physical meaning can be assigned to them (sec also 
Q)- Using equal time field correlation functions one can 
define a gluon number density n^kx), but the definition is 
not unique. The question of defining the number density 
is discussed in the following section. 



III. PARTICLES IN A CLASSICAL FIELD 



so we can identify: 



\n(k)\ 2 =u)(k)n(k), \<p(k)\ 2 = 



n(k) 
^(fc)' 



(25) 



For an interacting theory one can define the number 
distribution as follows: 



V m k )\ 

There is also another possibility, we can also assume a 
dispersion relation ujf rcc (k) = \f m 2 + k 2 and define 



n(k) 



KWl 2 ^ 

LO ilcc (k) ' 



(27) 



The latter approach is the one we take. Explicitly, for 
this particular theory described by the Hamiltonian, Eq. 
©, a 2-dimensional gauge field with an adjoint represen- 
tation scalar field on the lattice, we define 



x 2 1 



where 



9 -E?{k T )E?{-k T ) + ln"(k T )ir a (-k T ) 

(28) 



sin 



2 Qiky 



(29) 



is the free, massless lattice dispersion relation. We can 
then verify that our method is consistent with the ap- 
proach of Eq. (|26|l by looking at the correlation functions 



1 (E?(k T )E?(-k T )) 
rV (At(k T )Af(-k T )) 



and 



(7r°(fc T )7r (~fc T )) 



(r(k T )r(-k T )) 
( 3 °) 

and verifying that they behave as ui(k) ss k (see Fig. QJ. 



In a weakly interacting scalar theory it is easy to de- 
fine a particle number corresponding to a given classical 
field configuration. Take a free Hamiltonian and Fourier 
transform it: 



H = 



d d x 
d d k 



1 



<kt 



1 



2 (k)\m\ 2 



= / d d kuj(k)n(k), 



(22) 



(23) 



with the free dispersion relation ui 2 (k) = k 2 + m 2 . Av- 
eraged in time the energy is distributed equally between 
the degrees of freedom: 



\\-K(k)\* = \^{k)\4,{k)? 



(24) 



IV. RESULTS 

To state our results in a form easily comparable with 
Q let us define the same dimensionless quantities Jn and 
Je as follows: 



1 



dE, 



!e g^R\p? d v 



f. 



dN ir 



N 



g 2 irR 2 A p 2 drj 



(31) 
(32) 



As discussed in Sec. [n]the quantities $e and /at are func- 
tions of only one dimensionless variable g 4 irR A p 2 . In the 
weak field limit, namely for g^-nR 2 A p 2 < 50, Je and /at 
have a strong dependence on g A nR 2 A ^x 2 . This signals a 
dependence on the infrared cutoff of the theory. In the 
strong field limit, i.e. at large enough values of g 4 irR A p 2 , 
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FIG. 1: The functions IJ30J. The circles are uj(k) determined 
from the transverse fields E 1 and Ai, the solid line is uj(k) 
determined from ir and <f). The maximum value of ka is 2y/2. 



the nonlinearities of the infrared modes regulate this in- 
frared divergence and /e and /at become approximately 
independent of g 4 TrR A fj, 2 , as can be seen from Fig. 
Our results for the energy and multiplicity are summa- 
rized in Figs. Hand Hand Table D The total energy as a 
function of time in different field components is plotted 
in Fig. 0| and the energy in the different field components 
in Fig. 



FIG. 3: The functions f E and f N defined by Eqs. ||3TJ, J32J| 

for constant yj g 4 [i 2 nR 2 A and with different lattice spacings. 
The horizontal axis is g fia, so the continuum (a — * 0) limit is 
obtained by extrapolating each set of points to the g 2 (ia — 0- 
axis on the left. 





H (GeV) 


fE 




72 


0.29 


0.265 ± 0.005 


0.297 ± 0.006 


120 


0.49 


0.227 ± 0.003 


0.289 ± 0.003 


192 


0.78 


0.238 ± 0.005 


0.329 ± 0.006 



TABLE I: The values for /jv and fs corresponding to the 
points nearest to the continuum limit in Fig. The value 



0.7- < 

0.6- T 
_jz0.5- 
H 0.4- 1 

0.3- 

0.2 







50 



100 

/ 4 2 D 2 \' /2 

(8 li nR A ) 



150 



200 



FIG. 2: The functions /b and /at as defined by Eqs. 1)31^ . 
1321 1 vs. v/ g 4 H 2 ' 7 fR A . Computed on a 256 2 -lattice. 

Our result for Je is smaller than that of |4( by ap- 
proximately a factor of two. The function dN/ d 2 kT(kr) 
we obtain is different although its integral over fcr-space, 



148 fm . The value of /e is obtained by fitting the energy to 
a form A + Be~ T ^ T ° and using the value A. The multiplicity 
is measured at a time r = 5/fi, but its dependence on r is 
very weak. 



and thus /at, happens to be the same. This difference is 
illustrated in Fig. 

One can also derive a large kx analytic expression for 
the multiplicity in the classical field model Q (see also 
|10|). An expansion to the lowest nontrivial order in the 
field strength gives: 



diV 



1 N C {N C 



dry d 2 kx (2ir) 3 ir 



- l)gV fcr 2 
J A 2 



(33) 



with A some infrared cutoff. A useful check of the nu- 
merical computations is that they should approach the 
analytic result in the weak field limit of small g 2 fiR A , 
although the uncertainty from the infrared divergence 
of the analytical result can be numerically large. Fig- 
ure [7| shows that we do indeed observe a transition to 
a perturbative l/k^y. logarithmic factors - behaviour 
around kx > 2g 2 //, although in this region the shape 
of the spectrum is already severely modified by lattice 
effects, as can be seen comparing the plots for the two 
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FIG. 4: Total energy per unit rapidity as a function of time 
for n = 0.5 GeV. The three curves give an error estimate 
from 5 trajectories on a 512 2 -lattice. 
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FIG. 5: Energy in different field components from the same 
simulations as in Fig. 2] 



lattice sizes. But, as seen in Fig. [SJ the overall nor- 
malisation of our numerical result is far away from the 
analytical result at large \J g 4 fi 2 irR\ and approaches it 
only for y/g 4 fj, 2 TrR\ < 10. This would suggest that the 
weak field approximation used to obtain the analytical 
result ll-'i-il) is unsuitable for a quantitative understanding 
of this classical field model, whose justification lies, after 
all, in the argument of strong fields. 

Our definition of the gluon number spectrum, Eq. I|28|l , 
is based on equal time correlators of fields. These corre- 
lators are gauge dependent, which is a fundamental dif- 
ficulty in defining a multiplicity of gluons for this clas- 



as a function of k/g 2 fi for v/ <7 4 /i 2 7ri?A2 

>2 



120. The solid line is our result for a 512 -lattice, the dotted 
line for a 256 2 -lattice and the dashed line a fit to the numerical 
result of Q- The area under the curves (which is just /jv 
defined in Eq. jS3 ) is approximately the same (although the 
logarithmic scale makes this hard to see). The dashed curve 
practically falls on top of the solid one if the vertical axis is 
scaled by 2 and the horizontal by 1/2 — a signal of a difference 
in the normalisation. 
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FIG. 7: 



-JJC 

„<•,, -If 



dJV 



: r . , ^ as a function of k/g 2 fi from the same 

simulations as Fig. |S| The solid line is our result for a 512 2 - 
lattice and the dotted line for a 256 2 -lattice. 



sical field model. We have studied this gauge depen- 
dence by using as an example the electric field correlator 
{E?(k T )E?(-k T )), which is plotted in Fig. M Its gauge 
dependence is limited mainly by the constraint that the 
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FIG. 



1 dN 



as a function of ka plotted for 

^/g 4 fi 2 TTR 2 A = 240/i/GeV with values of given in the figure, 
compared with the analytical continuum result, Eq. 1331 . 



integral 



J A 2 k T Et{k T )Et{-k T ), 



(34) 



which is proportional to the energy in the electric field, 
is gauge independent. 

To determine the multiplicity using Eq. (|28[) we take 
the fields resulting from the initial conditions, Eqs. 114|) 
and evolve them in time according to the equations of 
motion, Eqs. ©. The "no gauge fixing" -curve in Fig. [5] 
shows the (Ef(kT)Ef (— fc^)) -correlator obtained in this 
way. The fields are then gauge transformed into the 
two dimensional Coulomb gauge diAi = to get the 
"Coulomb gauge" -correlator, also plotted in Fig. |5J This 
is the one that is used to determine the multiplicity. In 
Fig. [§] we also plot the same correlator in two other 
gauges, d x A x — and a "Coulomb + random" gauge, 
which is obtained by taking a field configuration in the 
Coulomb gauge and perfoming an independent random 
gauge transformation on each lattice site. In the latter 
the independent (Gaussian in this case) transformations 
on each lattice site naturally enhance the high momen- 
tum parts of the spectrum. 

According to the discussion in Sec [V] the value of the 
parameter /i relevant to RHIC phenomenology would 
be (x = 0.5 GeV or A s = 2 GeV. One can then ask 
whether this is indeed in the domain of validity of the 
present model, i.e. whether the occupation numbers of 
gluons are high enough. To address this question we 
plot in Fig. 1101 the two dimensional phase space density 

/(M = 2(^-1) ^"^7' where the s P in and color de- 
generacy has been divided out. It is of order one only up 
to momenta of a fraction of g 2 [i, meaning that the as- 



FIG. 9: The correlator {Et{k T )Ef{-k T )) in different gauges: 
the correlator resulting from the initial conditions and the 
equations of motion without additional gauge fixing, in "par- 
tial Coulomb" d x A x = gauge, in the Coulomb gauge 
diAi — and in a gauge obtained by a random gauge trans- 
formation of the Coulomb gauge field. 



sumption of high occupation numbers is only marginally 
satisfied. 

Seeing that the results of 0, are in many aspects 
qualitatively similar to ours and after a comparison of the 
numerical methods it seems that the difference in our re- 
sults concerning the energy and the number spectrum 
are simply due to a different normalisation of the SU(3) 
generators (compare Eq. © and Eq. (A5) of 11]). Any 
phenomenological discussion, such as the following, can- 
not be considered as an argument for the correctness of 
one or the other numerical result. 



V. PHENOMENOLOGY 

A. What to expect 

To discuss the phenomenological implications of these 
results in the light of RHIC experiments 0| one must 
relate the calculated initial multiplicities and transverse 
energies to the observed quantities. There are several 
scenarios that can be used to do this. Let us compare dif- 
ferent results with the assumption of early thermalisation 
and adiabatic expansion that has been successful in ex- 
plaining particle yields and elliptic flow. In this scenario 
the initial and final multiplicities, related by entropy con- 
servation, are approximately equal, and we take the total 
(charged and neutral) multiplicity per unit rapidity to be 



diVini 

dr/ 



dA fi 



rial 



dr] 



1000. 



(35) 
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conservation or some other mechanism would be 
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FIG. 10: The two dimensional phase space density f(k) 



(271-)-* dN 



2(JV C 



-1) KRi 



v as a function oik/g ji for ji = 0.5 GeV and 



g — 2, i.e. v/ 9 4 ^ 2 ^R 2 a — 120. The solid line is our result, the 
dashed line a fit to the numerical result of B. 



The observed transverse energy is 



d£ fi 



inal 



drj 



600 GeV. 



(36) 



The initial energy is larger than this, due to the expan- 
sion of the system. In a freely streaming system the en- 
ergy per unit rapidity is constant, whereas adiabatic lon- 
gitudinal expansion makes it decrease as t -1 / 3 . In ^| 
the energy is found to be reduced by a factor of 3.5. Be- 
cause the calculation of is done assuming a very early 
thermalization it gives an upper bound to the reduction. 
This translates into a bound for the initial transverse en- 
ergy d d'r"" ~ 2100 GeV. Thus a conservative estimate 
assuming "parton-hadron duality" , be it from entropy 



^<2.1GeV^. 
drj drj 



(37) 



The final state saturation model of is a pQCD- 
calculation supplemented by a sharp infrared cutoff de- 
termined from a simple geometrical final state saturation 
argument. The result of the calculation is (Eq. (5) of 



Psa t^± = 0.288 GeV 1 
drj 



1.050 



0.574 



Setting diVinit/dr/ to 1000 and taking A = 200, 
130 GeV gives p sat = 1.23 GeV. Then, from Eq. 
[l3| . we get 



(38) 
(7) of 



dE; 



init 



drj 



l-43p sat 
1.76 GeV 



dN 
drj 
dN 
drj 



(39) 



Intrinsically such an unphysically sharp infrared cutoff 
should produce too large an average energy per particle, 
because there are no gluons with pt < p sa ,t in the model. 
The constant coefficient in front of (|38|) is determined by 
the parton distributions and is not fitted to match the 
RHIC data. The result (|39|l could thus be regarded as a 
theoretical upper bound on the initial energy. 



B. Classical Yang-Mills result 

Let us take from Table [I] the result for /i w 0.5 GeV, 
which is the value of [i that gives approximately the right 
multiplicity. We get 



diV ir 



2 ,,2 



dr\ 



0.2V7Ti?^ 



This gives /i = 0.48 GeV or A s = 1.9 GeV. 
energy is 



dgini 
dr/ 



0.23 5 4 7ri?> 3 

dTVinit 



0.79 5 >- 



1.5 GeV 



d.77 

diVinit 

drj 



(40) 
Then the 

(41) 
(42) 
(43) 



This is well within the bound (|37() . 

The result of 4] is f N = 0.3. Setting diV init / d^ = 1000 
this gives us Agi?^ = 65. Taking ttR\ — 148 fm 2 this 
means A s = 1.87 GeV. For f E the result in is Je = 
0.537 for A S R A = 25 and f E = 0.497 for A S R A = 83.7. 
Taking the value Je = 0.5 one gets 



dE u 



drj 



1.67A S 



3.1 GeV 



dAf ir 



drj 

dN init 
drj 



(44) 
(45) 
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Thus the average energy per particle is 3.1 GeV, which is 
very hard to reconcile with the estimate Q37II and forces 
one to either give up the assumptions behind that es- 
timate or conclude that RHIC energies are not in the 
domain of validity of the classical field model. One can 
indeed argue, as in [4|, that gluon number increasing pro- 
cesses lower the average energy per particle in the subse- 
quent evolution of the system. 

VI. FINITE NUCLEI 



For k ;§> g 2 p this approaches the original McLerran- 
Venugopalan model, but for k <C g 2 p the fluctuations are 
cut off as (p a (kT)p h (kT)) ~ k 2 Our results for the multi- 
plicity and energy as a function of centrality are plotted 
in Fig. ^] All data points have been produced with the 
same number, 10, of configurations, the larger errors seen 
using the original Gaussian weight function are a signal 
of its strong dependence on few infrared modes. The dis- 
crepancy in the ratio E/N between our results and those 
of Ij, [111 remains also in the finite nucleus case. 



It is easy to naively generalise the model to finite 
nuclei. The Gaussian distribution of the random color 
charges is argued to arise from a sum of independent fluc- 
tuating charges. Thus it is the variance (p a (xt) p b (Ut)) 
that should be proportional to the thickness of the nu- 
cleus: 

(p a (x T )p\yr)) = g 2 p 2 S XT , yT S ab T(x T - x To ) (46) 



with T(xt) ~ \J R\ — xt 2 (or some more sophisticated 
thickness function). Note that the normalisation of p, is 
different from the square nucleus case, here we fix it by 
the condition 



]T (p a (x T )p b (y T ))=S ab g 2 p 2 nR 



(47) 



One can then proceed as previously. But the prob- 
lem one encounters is that the colour fields generated 
by the sources have long Coulomb tails outside the nu- 
clei. In two dimensions the initial colour fields (J3J decay 
only logarithmically away from the nuclei. Physically the 
colour fields should decay at distances ~ 1/Aqcd due to 
confinement physics not contained in this model. 

The approach of Q and also advocated by [l5| . 
is to directly address this question by imposing colour 
neutrality of the sources at a length scale of the order of 
a nucleon radius. But it is also possible that a proper 
inclusion of saturation effects in the probability distribu- 
tion of the initial colour sources might cure this problem. 
Saturation does, after all, suppress the very long wave- 
length modes responsible for the long tails. 

Exploring the full implications of the "JIMWLK" 
renormalization group equation for heavy ion collisions 
is out of the scope of this work, but in the spirit of, e.g., 
fl(i| we have tried substituting the correlation function 
(|46|l with the following procedure. We take random vari- 
ables f a (xr) distributed as: 



(f a (x T )f b (yT))=S XT , yT S ab T(x T ). 



(48) 



The original McLerran-Venugopalan model, equation 
would be obtained with the choice p a (xr) — gpf a (xT)- 
Now we Fourier transform and take 



p a (k 7 



k 2 



5AM 



0V 



;f a (k T ). 



(49) 
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FIG. 11: The functions j^ max /at (circles) and 
angles) vs. iVpart • iV max « 375 is iV part corresponding to 
impact parameter b = 0. The open symbols are results calcu- 
lated with the original Gaussian weight function and the filled 
symbols with the saturation ansatz l|490 . The conversion from 
b to A?p a rt from 13 ■ 



In [13 it is said that "our [using 'color neutral' initial 
conditions] results may be quantitatively similar to RG 
evolved predictions". This can also be seen in Fig. 1 of 
|ll| , where the neutrality condition originally imposed at 
the scale Aqcd has an effect up to the scale g 2 /i, leading 
to a modification of the Gaussian weight function that 
is very similar to ours. It might thus turn out that at 
RHIC energies it is not yet possible to distinguish effects 
from two physically very different phenomena, confine- 
ment and saturation. 



VII. CONCLUSIONS AND OUTLOOK 

We have applied the classical field approach to heavy 
ion collisions and calculated the energy and number den- 
sities and the spectra of the gluons produced in the initial 
stages of the collisions. We have also extended the model 
to finite nuclei and experimented with a crude saturation- 
inspired modification of the original model. The gauge 
dependence of equal time correlators of the fields which 



9 



makes it difficult to define a gluon number density has 
also been investigated. A more practical difficulty in 
the model is that the phase space density of particles at 
RHIC might not yet be large enough to justify its use, i.e. 
the saturation scale might not be large enough compared 
to Aqcd- For hard modes whose phase space density is 
small one does not even expect a classical field approach 
to work, and the transition to a pQCD regime should be 
understood better. 

Further things that need to be investigated within this 
approach include the incorporation of the "JIMWLK" 
renormalisation group equation into the calculation. A 
better understanding of thermalisation, if possible within 
the classical approach, might require extending the study 
to a 3+1-dimcnsional model. The "best estimate" in 



terms of physical postdictions for RHIC or predictions 
for LHC phenomenology is not settled yet, but is hope- 
fully converging. 
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